Continuum description of profile scaling in nanostructure decay 
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The relaxation of axisymmetric crystal surfaces with a single facet below the roughening transition 
is studied via a continuum approach that accounts for step energy gi and step-step interaction energy 
33 > 0. For diffusion-limited kinetics, free-boundary and boundary-layer theories are used for self- 
similar shapes close to the growing facet. For long times and ^ < 1, (a) a universal equation 

is derived for the shape profile, (b) the layer thickness varies as i^) 1 ^ 3 , ( c ) distinct solutions are 

found for different and (d) for conical shapes, the profile peak scales as (^■) _1/ ' 6 - These results 
compare favorably with kinetic simulations. 



o 

CZ3 



-a 
c 

o 
o 



> 
in 

00 

o 

O 
m 
o 

-I— > 



i 

-a 
c 

o 
o 



The drive toward smaller features in devices has fu- 
eled much interest in low-temperature kinetic processes 
such as growth, etching, and morphological relaxation. 
The constantly decreasing temperatures present increas- 
ing challenges for treatment of thermodynamics, kinetics 
and macroscopic evolution of surfaces. A crystal sur- 
face at equilibrium undergoes a roughening transition at 
a surface orientation-dependent temperature Tr 0, 0- 
In equilibrium at temperature T, crystal facets (planar 
regions of the surface) have Tr > T, whereas orienta- 
tions in continuously curved portions of the surface have 
Tr < T. In numerous nonequilibrium situations below 
Tr, a crystal surface relaxes to its equilibrium shape via 
the lateral motion of steps at a rate limited mainly by 
the diffusion of adatoms across terraces and attachment 
and detachment at step edges. Here we report a contin- 
uum treatment of this evolution using a partial differen- 
tial equation (PDE) and obtain scaling laws and universal 
aspects of the solutions. 

Morphological equilibration for surfaces above Tr is 
described by a continuum treatment 0, 0] where the sur- 
face free energy, which is an analytic function of the sur- 
face slope, and chemical potential Q are ingredients in a 
fourth-order PDE for the evolution of the surface profile. 
However, this analysis is not applicable to surfaces below 
Tr because the surface free energy is not analytic in the 
surface orientation j^ElS]; see Eq. @ below. 

Efforts to describe morphological evolution below Tr 
began in the mid 1980s and include simulations of the 
motion of monatomic crystalline steps and continuum 
thermodynamic approaches. In the latter, to account 
for evolution due to the motion of steps separating ter- 
races below the basal plane's Tr, the step density, which 
is proportional to the surface slope, is introduced as a 
variable within a coarse-grained continuum description 
on a scale large compared to the step separation (typ- 
ically, 1-10 nm). Expressions have been developed for 
the chemical potential of atoms at interacting step edges, 
leading to a nonlinear PDE for evolution of periodic sur- 



face modulations la, |9J; progress has been made toward 
solving the PDE |l0l . llll b ut has been hindered when evo- 
lution involves facets [Til Ilflj . Kinetic simulations mimic 
nanoscale processes and so have been used to describe the 
detailed motion of many steps, as well as the evolution 
of facets, via coupled differential equations [Til fljl llif. 
Nevertheless, the kinetic simulations are generally lim- 
ited in their ability to characterize universal features of 
the shape evolution. Here we show that the shape profile, 
including the facet, can be treated using a continuum, 
thermodynamic description that illuminates scaling as- 
pects of the kinetic behavior; for this purpose, we use an 
analytical framework that transcends the limitations of 
continuum approaches previously recognized [l^ . 

Our analytical approach treats facet evolution as a 
free-boundary problem 0, [^. The surface height is 
h(r,t), where r = (x,y) = re r is the position vector in 
the plane of reference; see Fig. 1 for an axisymmetric 
shape. The step density is |V/i|; V/i = (h x ,h y ) and sub- 
scripts denote partial derivatives. Denoting the atomic 
volume by fl and the surface current (atoms per length 
per time) by j, the conservation equation for adatoms is 



o 



(i) 



The current is j = — c s D s \7 /i/ksT, where D s and c s 
are the surface diffusivity and adatom concentration, and 
/i(r) is the chemical potential of atoms at step edges; D s 
is in principle a tensor function of V/i. 

We focus on diffusion-limited (DL) kinetics, where dif- 
fusion of adatoms across terraces is the rate-limiting pro- 
cess, and further assume that c s is constant and D s is a 
scalar constant. Equation becomes 

dh c SJ D s n vV (2) 



dt 



k B T 



Next, /j, and h are related via the surface free energy per 
unit projected area, G. A common expression for G of 
vicinal surfaces for T < Tr assumes that G depends on 
the step density according to 0, 0, Hf| 



G(Vh)=g Q + 9l \Vh\ + -g 3 \Vh\ s 



(3) 
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FIG. 1: View of an axisymmetric surface profile, on both 
the macroscale and the nanoscale where the atomic steps are 
evident. The evolution of surface morphology is caused by 
the motion of steps. 



The go term represents the surface free energy of the 
reference plane, 31 is the step energy (line tension), and 
33, which accounts for step-step interactions, includes en- 
tropic repulsions due to fluctuations at the step edges and 
pairwise energetic interactions between adjacent steps. 
All go, g\ and 33 are temperature dependent and we con- 
sider repulsive interactions between steps, 173 > 0. 

The surface chemical potential is derived from G by 



the relation H 0, E3 
(JL. JL) " 

^ dh x ' dh y j 
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where 



S(Vh) 
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(4) 



The surface evolution equation follows by combination 
of Eqs. lHJ-Q. We use cylindrical coordinates to de- 
scribe the relaxation of axisymmetric shapes h = h(r,t) 
that are smooth along the surface outside the facet and 



dh 
r dr , 



have negative slope, ^ < (Fig. 1). Since Vh = e 
it is convenient to define the dimensionless step density 
or surface slope F(r,t) — — The surface then evolves 
according to a fourth-order nonlinear PDE for F, 



dF 
~dt 



9i dr 



r or 



(5) 



c 3 _D s Q 2 gi 
k B T 



has dimensions 



The material parameter B 
(length) 4 /time. 

Equation JSJ is supplemented with the initial condition 
F(r, 0) = —H'(r), where H(r) = h(r, 0) is the initial sur- 
face profile with the properties H'{r) = for r < W 
(the initial facet radius) and H'(r) < for r > W . 
Also, there are four boundary conditions applied at the 
facet edge, r = w(t). In particular, the height h and 
the current j are continuous at the facet edge. The lat- 
ter condition, along with rj — > as r — > 00, ensures 
that the total mass is conserved. A consequence of Eq. 
(JSJ) and the initial conditions is that no other facets are 
formed. Another condition is slope continuity at the facet 
edge, F(w,t) = (i.e., local equilibrium l2ll 1221 and it 
is also consistent with kinetic simulations |14|1. It is 
shown below that Eq. (J5J furnishes F(r, t) = 0{\Jr — w) 



as r — > w + where the coefficient is time dependent [23j . 
Finally, for r < w, where there is a facet, we extend con- 
tinuously through Eqs. Q-@ the variable fj,, although 
it no longer represents the true chemical potential, and 



the quantity N = e r 7V : 



SG 
" S(Vh) 



whose divergence yields 



yit/fl 0, 01 • Hence the final two boundary conditions 
are continuity of /x and of N. Although we now have a 
mathematically complete set of boundary conditions for 
Eq. ©, the issue of the boundary conditions remains a 
topic of discussion [IE 113 • Nevertheless, we expect that 
the main analytical and scaling ideas given below are in- 
dependent of the detailed form of these conditions. 

The boundary conditions described above relate ^ to 
the derivatives of F 2 at r = w + , the facet height hf(t) 
and the facet radius w(t). By differentiating hf(t) — 
h(w + (t),t) in time, we deduce that, at r — w, 

B {\ - ^w[{F 2 )' - 2w(F 2 y - w 2 {F 2 )'"] \} = h f w 3 , (6) 
3i 

where the dot denotes the time derivative. Then, an 
examination of /i and N and their continuous extensions 
on the facet gives two more conditions at r = w [iij] : 



w[3{F 2 )' -w 2 {F 2 ) 
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w[3(F 2 )' - w(F 2 )"]. (7) 



We now treat Eq. © with conditions JHJ and (JJJ) and 

F = at the facet edge as a free-boundary problem ^3] : 
there is a moving facet for r < w(t), where F = 0, and 
this facet connects smoothly to the rest of the profile for 
r > w(t). Note that this problem statement is valid for 
arbitrary In general, there exist an "outer" region, 
where only the line-tension energy g\ is important, and 
an "inner" region in the neighborhood of the facet edge, 
where the step-step interaction energy 33 becomes signif- 
icant. Motivated by kinetic simulations with ^ < 1 |14j . 
we set e = |^ and treat Eq. J5j analytically for the case 
with small e, i.e., weak repulsive interactions. Because 
the small parameter e multiplies the highest-order spatial 
derivative in Eq. JSJ, the shape evolution can be treated 
with boundary-layer theory [2jJ. We start with the so- 
lution for e = where the corresponding facet radius 
w(t;e) is denoted w(t;0) — Wo(t). From Eq. JSJ), the 
zcroth-order outer solution F(r, t;0) = Fo(r,t) satisfies 
dFo/dt = 3B/r 4 , which is integrated subject to the ini- 
tial condition Fo(r, 0) = —H'(r) to give 



F (r,t) = ZBtr- 4 - H'(r) 



r > w (t) . (8) 



At the facet edge, F (wo 7 t) = iBt/wo — H'(wo) 7^ 0, so 
the slope profile is discontinuous; this feature motivates 
a singular perturbation analysis. 

The next step is to examine how the inclusion of a 
nonzero e renders the slope continuous by retaining the 
highest derivative in Eq. JSJ. We therefore consider a re- 
gion of width S(t; e) <C w in the neighborhood of the mov- 
ing facet edge, and describe the solution in this region in 
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FIG. 2: Numerical solutions of Eq. Hill with the 

boundary conditions /o(0) = and /o(oo) = 1. Curves a-e 
are parametrized by (ci,c 3 ) = (1.5, -.8183548), (2,-1.113031), 
(3,-1.72107502), (3.5, -2.0302102), (3.6, -2.09232155) and 
correspond to e = 9.2 x 10~ 3 ,1.9 x 10 -3 , 1.7 x 10 -4 , 
6.8 x 10~ 5 ,5.7 x 10" 5 . Inset: The dashed curves are de- 
scribed by Eq. 11311 for a conical initial shape and different e, 
while the solid curve shows C3 as a function of c\ from the 
numerical solutions of Eq. I1H . 

terms of the local variable rj = (r — w)/S. Thus, we seek 
a long-time similarity solution that depends on the dis- 
tance from the facet edge and time, F(r, t; e) = T{r\, t; e). 
We anticipate that, to leading order in e, 

F( V ,t) ~ a {t) f o{m e) , V =[r-w{t;e)]/S(t;e) , (9) 

where fo depends implicitly on e through the boundary 
conditions. Substitution of (0 into JSJ and balance of 
the leading-order terms in e gives 

Thus, must be time-independent and we take it 

equal to unity without affecting observable quantities 
such as F or w. It follows that S(t;e) = 0(e 1/3 ), inde- 
pendent of the (axisymmetric) initial conditions, which 
is a prediction for a scaling law for the boundary-layer 
width in the case of DL kinetics. The neglected terms 
in Eq. lfTo"|> are (^(e 1 / 3 ) <C 1. Furthermore, the leading- 
order facet radius is w (t) 4 ~ AB f*dt'w(t') 3 a (t), where 
w = f and A(t) = e-VStf^e). 

We next examine solutions of Eq. (|10|l along with the 
prescribed boundary conditions. First, this equation is 
integrated once via matching T(rj, t) with the outer solu- 
tion JSJl, i.e., taking 77 3> 1 and r — > w + simultaneously. 
We find a (t) = 3Btw- i ~H'(w) 1 which by Eq. © deter- 
mines the explicit time dependence of the surface slope. 
Because at this point we have imposed no restrictions 
other than axisymmetry on the initial shape, we have in 
fact obtained a universal equation for fo(r]), i-e., 

{fir = h - 1 , (11) 



which is to be solved with /o(0) = and fo{r] — > 00) = 1. 
Near the origin, fo(f]) has the behavior 

fo(v) ~ W/ 2 + c 3 ^ 3/2 + c 5V 5 ^ + c 6V 3 + (12) 

where all coefficients c n with n > 5 are known in terms 
of ci and C3. Equation (|llfl has a growing mode for 
•q 1, which must be suppressed in order to satisfy the 
condition at 00; thus, C3 is found numerically in terms of 
ci . We solve Eq. pi(l numerically and so obtain a family 
of similarity solutions fo (r)) for different values of c\ |2j| ; 
see Fig. 2. We next show how the solution /o(f?) depends 
on e, which requires imposing conditions such as J7J). 
The substitution of Eq. © into Eq. J7J and use of the 

relations (/o)^ =0 = c \ and (/o%' = 4c i c 3 from Ec l- 03 
yield two parametric equations for c\ and c 3 . In the case 
with a conical initial shape, discussed at length below, 
continuity of the variable \i implies [24| 

(cic 3 )e^ = ~ m4 + fD- 1 ] 173 . (13) 

The intersection of the curve (|13f) with the set of points 
(01,03) that result from numerically solving Eq. I|ll|) is 
shown in the inset of Fig. 2, and determines a value of 
e for each of the solution curves of the main part of the 
figure. Thus, we have determined a family of e-dependent 
similarity solutions fo (77; e) . 

There is one more scaling law that comes from the 
analysis. Each of the curves fo(r]) in Fig. 2 has a well- 
defined absolute maximum. Using l|12|) each maximum 
may be estimated to be 0(c^ 2 |c3| -1 / 2 ) and to occur at 
Vmax = 0{ci/\c 3 \), which is independent of e to leading 
order. Thus, according to Eq. (JT^J, C\ and c 3 are 0(e -1 / 6 ) 
and so the maximum slope is predicted to be 0(e -1 / 6 ). 

We now compare the predictions from this continuum 
approach based on Eq. (JSJ with the kinetic simulations 
for the DL case reported by Israeli and Kandel f° r a 
conical initial shape. In their simulations these authors 
vary a parameter, g = (const.) • 53, holding g\ fixed, 
which in our analysis is equivalent to changing e. Their 
simulations furnished a ^-dependent family of solutions 
(see their Figs. 4(b) and 6), which correspond to our e- 
dependent curves fo(i]\ e )- Israeli and Kandel also derived 
a complicated, (/-dependent differential equation in the 
scaling variable x — rf(Bt) 1 / 4 (not to be confused with 
the cartesian coordinate), which, in the limit of small g, 
effectively reduces to our (lllfl . However, on the basis of 
their equation they found multiple solutions whereas we 
provide a unique solution fo(i]; e) for each e. 

Next, we consider scaling behavior with e. First, we 
examine the scaling of the boundary layer near the facet 
edge. We define the boundary-layer thickness as the 
distance from the facet edge,a;o, to the position of the 
peak, Xpeak, of the step density F pce ^. In Fig. 3 we show 
the results of kinetic simulations (symbols) for a; pca k — xo 
vs. g and compare with our e 1 / 3 scaling prediction (solid 
line). Second, in Fig. 3 we examine how -F pG ak varies with 
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valuable feedback and for sharing with us details of their 
simulation results. 



step-step interaction strength, g 

FIG. 3: Log-log plot of the boundary-layer thickness 5(t; e) 
and the maximum of step density -Fp Ca k as functions of e. 
The crosses represent the results of kinetic simulations given 
to us by Israeli and Kandel |l4j| for the DL case. Here, 5(t; e) 
is estimated as the distance a: pea k — xo between the facet 
edge, xq = w{Bt)~ 1 ^ i , where F = 0, and the position :r pca k of 
the maximum of F. The straight lines correspond to the e 1 ' 3 
and e~ 1/l6 scaling laws predicted according to Eq. UUI . 



g, for which the results of kinetic simulations (symbols) 
are compared with the e -1 / 6 scaling prediction (solid 
line). In both cases the agreement is very good. With 
regard to the deviations in the boundary-layer width for 
g < ICR 6 , as e decreases in the simulations x pca k ap- 
proaches the facet edge so that the boundary-layer width 
is relatively small on the scale of the step spacing and 
is consequently poorly defined; its evaluation in discrete 
simulations thus becomes prone to errors. 

As shown above, qualitative predictions, such as the 
form of the multiple solution curves, and quantitative 
predictions, such as the e 1 / 3 scaling of the boundary- 
layer width and the e -1 / 6 scaling of the maximum of the 
slope, can be deduced from a continuum approach based 
on Eq. JSJ with the use of free-boundary and boundary- 
layer theories |29| . Further, simple analytical arguments 
show that, for any admissible initial slope F(r, 0) = fix", 
for a wide range of v including — 4 < v < 1, the facet 
radius is w = 0{t 1 ^ l/+A ^) at sufficiently long times. In 
addition, we expect that, for a class of non-axisymmetric 
initial shapes, the near-facet boundary layer width still 
retains the 0(e 1//3 ) scaling for the isotropic surface free 
energy of Eq. J3J). 

The continuum approach and the free-boundary view- 
point capture the essential physics of crystal surface evo- 
lution below Tr. This development should give further 
impetus to continuum approaches to morphological evo- 
lution even at the nanoscale for structures far below Tr. 
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